#include<math.h>
#include<iostream>
using namespace std;
double abs1(double n)
{
    if(n<0)
        n=-n;
    return n;
}
double *zhuanzhi(double *a,int m,int n)//矩阵转置 a为m行n列
{
    double *c=new double[m*n];
    for(int i=0;i<m;i++)
    {
        for(int j=0;j<n;j++)
        {
            c[m*j+i]=a[n*i+j];
        }
    }
    return c;
}
double *operatejia(double *a,int m,int n,double *b)//矩阵加法 a,b均为m行n列
{
    double *c=new double[m*n];
    for(int i=0;i<m*n;i++)
    {
        c[i]=a[i]+b[i];
    }
    return c; 
}  
double *operatejian(double *a,int m,int n,double *b)//矩阵减法 a,b均为m行n列
{
    double *c=new double[m*n];
    for(int i=0;i<m*n;i++)
    {
        c[i]=a[i]-b[i];
    }
    return c; 
}
double *operateshucheng(double *a,int m,int n,double x)//矩阵数乘 a为m行n列,数乘为x
{
    double *c=new double[m*n];
    for(int i=0;i<m*n;i++)
    {
        c[i]=a[i]*x;
    }
    return c; 
}
double *operatechengfa(double *a,int m,int n,int x,double *b)//矩阵乘法 a为m行n列,b为n行x列
{
    double *c=new double[m*x];
    for(int i=0;i<m;i++)
    {
        for(int j=0;j<x;j++)
        {
            c[x*i+j]=0;
            for(int k=0;k<n;k++)
            {
                c[x*i+j]+=a[n*i+k]*b[x*k+j];
            }
        }
    }
    return c; 
}

void lowmatrix(double a[], int n,double b[])//前代法 针对方阵
{ 
    for(int j=0;j<n-1;j++)
    {
        b[j]=1.0*b[j]/a[j+n*j];
        for(int i=j+1;i<=n-1;i++)
        {
            b[i]=b[i]-b[j]*a[j+n*i];
        }
    }
    b[n-1]=b[n-1]/a[n*n-1];
}
void uppermatrix(double a[], int n,double b[])//回代法 针对方阵
{

    for(int j=n-1;j>0;j--)
    {
        b[j]=b[j]/a[j+n*j];
        for(int i=0;i<=j-1;i++)
        {
            b[i]=b[i]-b[j]*a[j+n*i];
        }
    }
    b[0]=1.0*b[0]/a[0];
}
void LU(double a[],int n)//LU分解 针对方阵
{
    for(int i=0;i<n-1;i++)
    {
        for(int j=i+1;j<n;j++)
        {
            a[n*j+i]=a[n*j+i]/a[n*i+i];
            for(int k=i+1;k<n;k++)
                {
                    a[n*j+k]=a[n*j+k]-a[n*j+i]*a[n*i+k];
                }
        }
    }              
}
void lieLU(double a[],int n,double b[])//列主元LU分解 针对方阵
{
    for(int i=0;i<n-1;i++)
    {
		int p=i;
		for(int k=i+1;k<n;k++)
		{
	    	if(a[n*k+i]*a[n*k+i]>a[n*p+i]*a[n*p+i])
	    	{
				p=k;
	    	}
    	}
        for(int j=0;j<n;j++)
        {
	    	double c=a[n*p+j];
	    	a[n*p+j]=a[n*i+j];
	    	a[n*i+j]=c;
		}
        for(int j=0;j<n;j++)
		{
            double c=b[n*p+j];
            b[n*p+j]=b[n*i+j];
            b[n*i+j]=c;
  		}
		if(a[n*i+i]!=0)
		{
	    	for(int j=i+1;j<n;j++)
            {
                a[n*j+i]=a[n*j+i]/a[n*i+i];
                for(int k=i+1;k<n;k++)
                {
                    a[n*j+k]=a[n*j+k]-a[n*j+i]*a[n*i+k];
                }
	    	}
		}
    }
}
void cholesky(double a[],int n) //平方根法
{
    for(int i=0;i<n;i++)
    {
		a[n*i+i]=sqrt(a[n*i+i]);
		for(int p=i+1;p<n;p++)
		{
	    a[n*p+i]=a[n*p+i]/a[n*i+i];
		}
		for(int j=i+1;j<n;j++)
		{
	    	for(int k=j;k<n;k++)
	        	a[n*k+j]=a[n*k+j]-a[n*k+i]*a[n*j+i];
		}
    }
    for(int i=0;i<n;i++)
    {
        for(int j=i+1;j<n;j++)
        {
	    	a[n*i+j]=0;
        }
    }
}
void gaijincholesky(double a[],int n,double b[]) //改进的平方根法
{
    for(int i=0;i<n;i++)
    {
		for(int j=0;j<i;j++)
        {
	    	b[j]=a[n*i+j]*a[n*j+j];
		}
        for(int j=0;j<i;j++)
		{
	    	a[n*i+i]=a[n*i+i]-a[n*i+j]*b[j];
	    	for(int k= i+1;k<n;k++)
	    	{
				a[n*k+i]=a[n*k+i]-a[n*k+j]*b[j];
	    	}
		}
	    for(int k=i+1;k<n;k++)
			a[n*k+i]=a[n*k+i]/a[n*i+i];
    } 
    for(int i=0;i<n;i++)
    {
        for(int j=i+1;j<n;j++)
        {
	    a[n*i+j]=0;
        }
    }
}
double *qiuni(double *a,int n)
{
    double *d=new double[n*n];
    double c[n];
    LU(a,n);
    for(int i=0;i<n;i++)
    {
		for(int j=0;j<n;j++)
		{
	    	if(i==j)
				d[n*i+j]=1;
	    	else if(i>j)
				d[n*i+j]=a[n*i+j];
	    	else
				d[n*i+j]=0;
		}
    }
    for(int i=0;i<n;i++)
    {
		for(int j=0;j<n;j++)
		{
	    	if(j==i)
	        	c[j]=1.0;
	    	else
				c[j]=0;
		}
        lowmatrix(d,n,c);
        uppermatrix(a,n,c);
        for(int j=0;j<n;j++)
            d[n*j+i]=c[j];
    }	
    return d;
} 
void jiefangcheng(double a[],int n,double b[])
{
    double c[n*n];
    for(int i=0;i<n*n;i++)
    {
		c[i]=a[i];
    }
    double d[n*n];
    for(int i=0;i<n;i++)
    {
		for(int j=0;j<n;j++)
		{
	    	if(j==i)
	        	d[n*i+j]=1;
	    	else
			d[n*i+j]=0;
		}
    }
    lieLU(c,n,d);
    double *e=new double[n*n];
    e=operatechengfa(d,n,n,1,b);
    double f[n*n];
    for(int i=0;i<n;i++)
    {
		for(int j=0;j<n;j++)
		{
	    	if(i==j)
				f[n*i+j]=1;
	    	else if(i>j)
				f[n*i+j]=a[n*i+j];
	    	else
				f[n*i+j]=0;
		}
    }
    lowmatrix(f,n,e);
    uppermatrix(c,n,e);
    for(int i=0;i<n;i++)
    {
		b[i]=e[i];
    }
    delete []e;
}
double sign(double n)
{   
    double s=0;
    if(n>0)
		s=1;
    if(n<0)
		s=-1;
    return s;
}
double *sign(double a[],int n)
{
    double *c=new double[n];
    for(int i=0;i<n;i++)
    {
		c[i]=sign(a[i]);
    }
    return c;
}
double xiangliangwuqiongfanshu(double a[],int n)
{
    double s=abs1(a[0]);
    for(int i=1;i<n;i++)
    {
        if(s<abs1(a[i]))
	     	s=abs1(a[i]);
    }
    return s;
}
double xiangliangyifanshu(double a[],int n)
{
    double s=0;
    for(int i=0;i<n;i++)
    {
		s+=abs1(a[i]);
    }
    return s;
}
double juzhengwuqiongfanshu(double a[],int n)
{
    double s=0;
    for(int i=0;i<n;i++)
    {
		s+=abs1(a[i]);
    }
    for(int i=0;i<n;i++)
    {
		int p=0;
		for(int j=0;j<n;j++)
        {
	    	p+=abs1(a[n*i+j]);
		}
		if(p>s)
	    	s=p;
    }
    return s;
}	
double juzhengnidezhuanzhideyifanshu(double a[],int n)
{
    double x[n];
    for(int i=0;i<n;i++)
        x[i]=1.0/n;
    int k=1;
    int p=0;
    double *f=new double[n*n];
    f=zhuanzhi(a,n,n);
    double v;
    while(k==1)
    { 
        double*b=new double[n];
		for(int i=0;i<n;i++)
		{
	    	b[i]=x[i];
		}
		jiefangcheng(f,n,b);
		double *e=new double[n];
		e=sign(b,n);
		double *z=new double[n];
		for(int i=0;i<n;i++)
		{
	    	z[i]=e[i];
		}
		jiefangcheng(a,n,z);
        double *d=new double[1];
		d=operatechengfa(z,1,n,1,x);
        if(xiangliangwuqiongfanshu(z,n)<=d[0])
		{
	    	v=xiangliangyifanshu(b,n);
	    	k=0;
	    	break;
		}
		else
		{
	    	for(int i=0;i<n;i++)
	    	{
	    	    if(xiangliangwuqiongfanshu(z,n)==abs1(z[i]))
				{
		 			p=i;
		    		break;
				}
	    	}
	    	for(int i=0;i<n;i++)
	    	{
	        	if(i==p)
		    		x[i]=1;
		 		else
		    		x[i]=0;
	    	}
		}
		delete []b;
		delete []e;
		delete []z;
		delete []d;
    }
    return v;
}
double *Householder(double a[],int n)
{
    double *x=new double[n];
    double p=0;
    double b=xiangliangwuqiongfanshu(a,n);
    double aa[n];
    for(int i=0;i<n;i++)
    {
		aa[i]=a[i];
    }
    for(int i=0;i<n;i++)
    {
		aa[i]=1.0*aa[i]/b;
    }
    double d=0;
    for(int i=1;i<n;i++)
    {
		d+=aa[i]*aa[i];
    }
    double c=sqrt(aa[0]*aa[0]+d);
    for(int i=1;i<n;i++)
    {
		x[i]=aa[i];
    }
    if(d==0)
    {
		x[0]=1;
        p=0;
    }
    else
    {
		if(aa[0]<=0)
	    	x[0]=aa[0]-c;
		else
	    	x[0]=(-d)*1.0/(aa[0]+c);
		p=2*1.0/((d+x[0]*x[0])*b*b);
		for(int i=0;i<n;i++)
		{
		    x[i]=x[i]*b;
        }
    } 
    return x;
}
double Householder1(double a[],int n)
{
    double x[n];
    double p=0;
    double b=xiangliangwuqiongfanshu(a,n);
    double aa[n];
    for(int i=0;i<n;i++)
    {
		aa[i]=a[i];
    }
    for(int i=0;i<n;i++)
    {
		aa[i]=aa[i]/b;
    }
    double d=0;
    for(int i=1;i<n;i++)
    {
		d+=aa[i]*aa[i];
    }
    for(int i=1;i<n;i++)
    {
		x[i]=aa[i];
    }
    double c=sqrt(aa[0]*aa[0]+d);
    if(d==0)
    {
		x[0]=1;
        p=0;
    }
    else
    {
		if(a[0]<=0)
		    x[0]=aa[0]-c;
		else
		    x[0]=(-d)/(aa[0]+c);
		p=2*1.0/((d+x[0]*x[0])*b*b);
		for(int i=0;i<n;i++)
		{
	    	x[i]=x[i]*b;
        }
    } 
    return p;
}
double *Householderjuzheng(double a[],int n)
{
	double aa[n];
	for(int i=0;i<n;i++)
	{
		aa[i]=a[i];
	}
	double *b=new double[n];
	b=Householder(aa,n);
	double c=Householder1(aa,n);
	double *d=new double[n*n];
	d=operatechengfa(b,n,1,n,b);
	double *e=new double[n*n];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			if(i==j)
				e[n*i+j]=1-c*d[n*i+j];
			else
				e[n*i+j]=-c*d[n*i+j];
		}
	}
	return e;
}
void QR(double a[],int m,int n,double b[])
{
	double d;
    for(int i=0;i<n;i++)
    {
		if(i<m-1)
		{
	    	double c[m-i];
	    	for(int j=0;j<m-i;j++)
	    	{
				c[j]=1.0*a[n*(i+j)+i];
	    	}
	    	double *x=new double[m-i];
	    	x=Householder(c,m-i);
	    	d=Householder1(c,m-i);
	    	double *f=new double[(m-i)*(m-i)];
	    	f=operatechengfa(x,m-i,1,m-i,x);
	    	double *e=new double[(m-i)*(m-i)];
	    	e=operateshucheng(f,m-i,m-i,d);
	    	double h[(m-i)*(m-i)];
	    	for(int j=0;j<m-i;j++)
	    	{
				for(int k=0;k<m-i;k++)
				{
		    		if(j==k)
						h[(m-i)*j+k]=1.0;
		    		else
						h[(m-i)*j+k]=0;
				}
	    	}
	    	double *g=new double[(m-i)*(m-i)];
	    	g=operatejian(h,m-i,m-i,e);
	    	double l[(m-i)*(n-i)];
	    	for(int j=0;j<m-i;j++)
	    	{
				for(int k=0;k<n-i;k++)
				{
		    		l[(n-i)*j+k]=a[(i+j)*n+i+k];
				}
	    	}
	    	double *o=new double[(m-i)*(n-i)];
	    	o=operatechengfa(g,m-i,m-i,n-i,l);
	    	for(int j=0;j<m-i;j++)
	    	{
				for(int k=0;k<n-i;k++)
				{
		    		l[(n-i)*j+k]=a[(i+j)*n+i+k];
				}
	    	}
	    	for(int j=0;j<m-i;j++)
	    	{
				for(int k=0;k<n-i;k++)
				{
		    		a[(i+j)*n+i+k]=o[(n-i)*j+k];
				}
	    	}
	    	for(int j=0;j<m-i;j++)
	    	{

		    	a[n*(i+j+1)+i]=x[j+1];
	    	} 
	    	b[i]=d;
	    	delete []x;
	    	delete []e;
	    	delete []f;
	    	delete []g;
	    	delete []o;
		}
    }
}
double *QR_R(double a[],int m,int n)
{
    double *x=new double[n*n];
    double b[n];
    for(int i=0;i<n;i++)
    {
		b[i]=0;
    }
    double aa[m*n];
    for(int i=0;i<m*n;i++)
    {
    	aa[i]=a[i];
    }
    QR(aa,m,n,b);
    for(int i=0;i<n;i++)
    {
		for(int j=0;j<n;j++)
		{
	    	if(i<=j)
				x[n*i+j]=aa[n*i+j];
	    	else
				x[n*i+j]=0;
		}
    }
    return x;
}
double *QR_Q(double a[],int m,int n)
{
	double aa[m*n];
	for(int i=0;i<m*n;i++)
	{
		aa[i]=a[i];
	}
    double *b=new double[m*m];
    double bb[m];
	for(int i=0;i<m;i++)
	{
		bb[i]=aa[i*n];
	}
	b=Householderjuzheng(bb,m);
	double *cc=new double[m*n];
	cc=operatechengfa(b,m,m,n,aa);
	for(int i=0;i<m*n;i++)
	{
		aa[i]=cc[i];
	}
    for(int i=1;i<n;i++)
    {
		double c[i*i];
		for(int j=0;j<i;j++)
		{
			for(int k=0;k<i;k++)
				{
					if(j==k)
						c[i*j+k]=1;
					else
						c[i*j+k]=0;
				}
		}
		double e[m-i];
		for(int j=0;j<m-i;j++)
		{
			e[j]=aa[n*(i+j)+i];
		}
		double *f=new double[(m-i)*(m-i)];
		f=Householderjuzheng(e,m-i);
		double h[m*m];
		for(int j=0;j<m;j++)
		{
			for(int k=0;k<m;k++)
			{
				if(j<i && k<i)
					h[m*j+k]=c[i*j+k];
				else if(j>=i && k>=i)
					h[m*j+k]=f[(m-i)*(j-i)+(k-i)];
				else 
					h[m*j+k]=0;
			}
		}
		double *o=new double[m*m];
		o=operatechengfa(b,m,m,m,h);
		for(int j=0;j<m;j++)
		{
			for(int k=0;k<m;k++)
			{
				b[m*j+k]=o[m*j+k];
			}
		}
		double *g=new double[m*n];
		g=operatechengfa(h,m,m,n,aa);
		for(int k=0;k<m*n;k++)
		{
			aa[k]=g[k];
		}
		delete []g;
		delete []o;
		delete []f;
	}
	return b;
} 
double* QRjiefangcheng(double a[],int m,int n,double b[])
{
	double *c=new double[m*m];
	c=QR_Q(a,m,n);
	double *d=new double[n*n];
	d=QR_R(a,m,n);	
	double *e=new double[m*n];
	for(int i=0;i<m;i++)
	{
		for(int j=0;j<n;j++)
		{	
			e[n*i+j]=c[m*i+j];
		}
	}
	double *f=new double[n*m];
	f=zhuanzhi(e,m,n);
	double *g=new double[n];
	g=operatechengfa(f,n,m,1,b);
	uppermatrix(d,n,g);
	return g;
}
double* jacobi(double a[],int n,double b[])
{
	double*c=new double[n];	
	for(int i=0;i<n;i++)
	{
		c[i]=1.0;
	}
    double x[n];	
	for(int i=0;i<n;i++)
	{
		x[i]=0;
	}
	for(int k=0;k<10000;k++)
	{
		for(int i=0;i<n;i++)
		{
			for(int j=0;j<n;j++)
			{
				if(i!=j)
					x[i]+=-1.0/a[n*i+i]*a[n*i+j]*c[j];
			}
			x[i]+=1.0*b[i]/a[n*i+i];
		}
		for(int i=0;i<n;i++)
		{
			c[i]=x[i];
		}
		for(int i=0;i<n;i++)
		{
			x[i]=0;
		}
	}
	return c;
}
double* GS(double a[],int n,double b[])
{
	double*c=new double[n];	
	for(int i=0;i<n;i++)
	{
		c[i]=1.0/n;
	}
    double x[n];	
	for(int i=0;i<n;i++)
	{
		x[i]=0;
	}
	for(int k=0;k<10000;k++)
	{
		for(int i=0;i<n;i++)
		{
			for(int j=0;j<n;j++)
			{
				if(i<j)
					x[i]+=-1.0/a[n*i+i]*a[n*i+j]*c[j];
				if(i>j)
					x[i]+=-1.0/a[n*i+i]*a[n*i+j]*x[j];
			}
			x[i]+=1.0*b[i]/a[n*i+i];
		}
		for(int i=0;i<n;i++)
		{
			c[i]=x[i];
		}
		for(int i=0;i<n;i++)
		{
			x[i]=0;
		}
	}
	return c;
}
double* SOR(double a[],int n,double b[],double w)
{
	double*c=new double[n];	
	for(int i=0;i<n;i++)
	{
		c[i]=1.0/n;
	}
    double x[n];	
	for(int i=0;i<n;i++)
	{
		x[i]=0;
	}
	for(int k=0;k<100;k++)
	{
		for(int i=0;i<n;i++)
		{
			x[i]+=(1-w)*c[i];
			for(int j=0;j<n;j++)
			{
				if(i<j)
					x[i]+=-w*1.0/a[n*i+i]*a[n*i+j]*c[j];
				if(i>j)
					x[i]+=-w*1.0/a[n*i+i]*a[n*i+j]*x[j];
			}
			x[i]+=w*1.0*b[i]/a[n*i+i];
		}
		for(int i=0;i<n;i++)
		{
			c[i]=x[i];
		}
		for(int i=0;i<n;i++)
		{
			x[i]=0;
		}
	}
	return c;
}
double xiangliangerfanshu(double a[],int n)
{
	double x=0;
	for(int i=0;i<n;i++)
	{
		x+=a[i]*a[i];
	}
	x=sqrt(x);
	return x;
}
double *gongetidu(double a[],int n,double b[])
{
	double *x=new double[n];
	for(int i=0;i<n;i++)
	{
		x[i]=1;
	}
	int k=0;
	double *c=new double[n];
	c=operatechengfa(a,n,n,1,x);
	double *d=new double[n];
	d=operatejian(b,n,1,c);
	double *e=new double[1];
	e=operatechengfa(d,1,n,1,d);
	double *p=new double[n];
	double f;
	double m;
	while(sqrt(e[0])>0.00000001*xiangliangerfanshu(b,n) && k<100000)
	{
		k=k+1;
		if(k=1)
		{
			for(int i=0;i<n;i++)
			{
				p[i]=d[i];
			}
		}
		else
		{
			f=e[0]/m;
			for(int i=0;i<n;i++)
			{
				p[i]=d[i]+f*p[i];
			}
		}
		double *g=new double[n];
		g=operatechengfa(a,n,n,1,p);
		double *h=new double[1];
		h=operatechengfa(p,1,n,1,g);
		double l=1.0*e[0]/h[0];
		for(int i=0;i<n;i++)
		{
			x[i]=x[i]+l*p[i];
		}
		for(int i=0;i<n;i++)
		{
			d[i]=d[i]-l*g[i];
		}
		m=e[0];
		e[0]=xiangliangerfanshu(d,n)*xiangliangerfanshu(d,n);
		delete []g;
		delete []h;
	}
	return x;
}	
double *juzhengduoxiangshi(double a[],int n)
{
    double *b=new double[n*n];
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            if(i==0)
               b[n*i+j]=-a[j];
            else if(i-j==1)
            	b[n*i+j]=1;
            else
            	b[n*i+j]=0;
         }
     }
	return b;
}
double mifa(double a[],int n)
{
	double *x=new double[n];
	x[0]=1;
	x[1]=1;
	for(int i=2;i<n;i++)
	{
		x[i]=0;
	}
	double c=0;
	for(int k=0;k<100000;k++)
	{
		for(int i=0;i<n;i++)
		{
			c=x[0];
			if(abs1(c)<abs1(x[i]))
				c=x[i];
		}
		for(int i=0;i<n;i++)
		{
			x[i]=x[i]/c;
		}
		double *d=new double[n];
		d=operatechengfa(a,n,n,1,x);
		for(int i=0;i<n;i++)
		{
			x[i]=d[i];
		}
	}
	double e=x[0];
	for(int i=0;i<n;i++)
	{
		if(abs1(e)<abs1(x[i]))
			e=x[i];
	}
	return e;
}
double *Hessenberg_Q(double a[],int n)
{
	double aa[n*n];
	double *q=new double[n*n];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			if(i==j)
				q[n*i+j]=1;
			else
				q[n*i+j]=0;
		}
	}
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			aa[n*i+j]=a[n*i+j];
		}
	}
	for(int k=1;k<n-1;k++)
	{
		double b[n-k];
		for(int i=0;i<n-k;i++)
		{
			b[i]=aa[n*(i+k)+k-1];
		}
		double *c=new double[(n-k)*(n-k)];
		c=Householderjuzheng(b,n-k);
		double d[n*n];
		for(int i=0;i<n;i++)
		{
			for(int j=0;j<n;j++)
			{
				if(i<k && j<k && i==j)
					d[n*i+j]=1;
				else if(i>=k &&j>=k)
					d[n*i+j]=c[(n-k)*(i-k)+j-k];
				else
					d[n*i+j]=0;
			}
		}
		double *e=new double[n*n];
		e=operatechengfa(q,n,n,n,d);
		for(int i=0;i<n;i++)
		{
			for(int j=0;j<n;j++)
			{
				q[n*i+j]=e[n*i+j];
			}
		}
		delete []c;
		delete []e;
	}			
	return q;
}
double *Hessenberg(double a[],int n)
{
	double aa[n*n];
	double *q=new double[n*n];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			if(i==j)
				q[n*i+j]=1;
			else
				q[n*i+j]=0;
		}
	}
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			aa[n*i+j]=a[n*i+j];
		}
	}
	for(int k=1;k<n-1;k++)
	{
		double b[n-k];
		for(int i=0;i<n-k;i++)
		{
			b[i]=aa[n*(i+k)+k-1];
		}
		double *c=new double[(n-k)*(n-k)];
		c=Householderjuzheng(b,n-k);
		double d[n*n];
		for(int i=0;i<n;i++)
		{
			for(int j=0;j<n;j++)
			{
				if(i<k && j<k && i==j)
					d[n*i+j]=1;
				else if(i>=k &&j>=k)
					d[n*i+j]=c[(n-k)*(i-k)+j-k];
				else
					d[n*i+j]=0;
			}
		}
		double *e=new double[n*n];
		e=operatechengfa(q,n,n,n,d);
		for(int i=0;i<n;i++)
		{
			for(int j=0;j<n;j++)
			{
				q[n*i+j]=e[n*i+j];
			}
		}
		delete []c;
		delete []e;
	}
	double *f=new double[n*n];
	f=operatechengfa(aa,n,n,n,q);
	double *g=new double[n*n];
	g=zhuanzhi(q,n,n);
	double *h=new double[n*n];
	h=operatechengfa(g,n,n,n,f);
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			if(i>=j+2)
				h[n*i+j]=0;
		}
	}
	return h;
}
double *shuangchongbuQR(double hh[],int n)
{
	int m=n-1;
	double *h=new double[n*n];
	for(int i=0;i<n*n;i++)
	{
		h[i]=hh[i];
	}
	double s,t,x,y,z;
	if(n==4)
	{
	s=h[n*(m-1)*(m-1)]+h[n*2-1];
	t=h[n*(m-1)*(m-1)]*h[n*2-1]-h[n*(m-1)+n-1]*h[n*(n-1)+m-1];
	x=h[0]*h[0]+h[1]*h[n]-s*h[0]+t;
	y=h[n]*(h[0]-s+h[n+1]);
	z=h[n]*h[2*n+1];
	}
	else
	{
	s=h[n*(m-1)+(m-1)]+h[n*2-1];
	t=h[n*(m-1)+(m-1)]+h[n*2-1]-h[n*(m-1)+n-1]*h[n*(n-1)+m-1];
	x=h[0]*h[0]+h[1]*h[n]-s*h[0]+t;
	y=h[n]*(h[0]-s+h[n+1]);
	z=h[n]*h[2*n+1];
	}
	for(int k=0;k<n-2;k++)
	{
		double xiangliang[3]={x,y,z};
		double *h0=new double[9];
		h0=Householderjuzheng(xiangliang,3);
		int q=max(1,k);
		double *h1=new double[3*(n-q+1)];
		for(int i=0;i<3;i++)
		{
			for(int j=0;j<n-q+1;j++)
			{
				h1[(n-q+1)*i+j]=h[n*(k+i)+q-1+j];
			}
		}
		double *h2=new double[3*(n-q+1)];
		h2=operatechengfa(h0,3,3,n-q+1,h1);
		for(int i=k;i<k+3;i++)
		{
			for(int j=q-1;j<n;j++)
			{
				h[n*i+j]=h2[(n-q+1)*(i-k)+j-q+1];
			}
		}
		int r=min(k+4,n);
		double h3[3*r];
		for(int i=0;i<r;i++)
		{
			for(int j=k;j<k+3;j++)
			{
				h3[3*i+j-k]=h[n*i+j];
			}
		}
		double *h4=new double[3*r];
		h4=operatechengfa(h3,r,3,3,h0);
		for(int i=0;i<r;i++)
		{
			for(int j=k;j<k+3;j++)
			{
				h[n*i+j]=h4[3*i+j-k];
			}
		}
		x=h[n*(k+1)+k];
		y=h[n*(k+2)+k];
		if(k<n-3)
			z=h[n*(k+3)+k];
			delete []h0;
			delete []h1;
			delete []h2;
			delete []h4;
	}
	double xiangliang2[2]={x,y};
	double *h5=new double[4];
	h5=Householderjuzheng(xiangliang2,2);
	double *h7=new double[6];
	for(int i=0;i<2;i++)
	{
		for(int j=0;j<3;j++)
		{
			h7[3*i+j]=h[n*(i+n-2)+n-3+j];
		}
	}
	double *h6=new double[6];
	h6=operatechengfa(h5,2,2,3,h7);
	for(int i=0;i<2;i++)
	{
		for(int j=0;j<3;j++)
		{
			h[n*(n-2+i)+n-3+j]=h6[3*i+j];
		}
	}
	double *h8=new double[2*n];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<2;j++)
		{
			h8[2*i+j]=h[n*i+j+n-2];
		}
	}
	double *h9=new double[2*n];
	h9=operatechengfa(h8,n,2,2,h5);
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<2;j++)
		{
			h[n*i+j+n-2]=h9[2*i+j];
		}
	}
	return h;
}
double *shuangchongbuQR_Q(double hh[],int n)
{
	int m=n-1;
	double *Q=new double[n*n];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			if(i==j)
				Q[n*i+j]=1;
			else
				Q[n*i+j]=0;
		}
	}
	double *h=new double[n*n];
	for(int i=0;i<n*n;i++)
	{
		h[i]=hh[i];
	}
	double s=h[n*(m-1)*(m-1)]+h[n*2-1];
	double t=h[n*(m-1)*(m-1)]*h[n*2-1]-h[n*(m-1)+n-1]*h[n*(n-1)+m-1];
	double x=h[0]*h[0]+h[1]*h[n]-s*h[0]+t;
	double y=h[n]*(h[0]-s+h[n+1]);
	double z=h[n]*h[2*n+1];
	for(int k=0;k<n-2;k++)
	{
		double xiangliang[3]={x,y,z};
		double *h0=new double[9];
		h0=Householderjuzheng(xiangliang,3);
		double *h00=new double[n*n];
		if(k==0)
		{
			for(int i=0;i<n;i++)
			{
				for(int j=0;j<n;j++)
				{
					if(i<3 &&j<3)
						h00[n*i+j]=h0[3*i+j];
					else if(i>=3 &&j>=3 &&i==j)
						h00[n*i+j]=1;
					else
						h00[n*i+j]=0;
				}
			}
		}
		else
		{
			for(int i=0;i<n;i++)
			{
				for(int j=0;j<n;j++)
				{
					if(i<3+k &&i>=k &&j<3+k && j>=k)
						h00[n*i+j]=h0[3*(i-k)+j-k];
					else if(i>=3+k &&j>=3+k &&i==j)
						h00[n*i+j]=1;
					else if(i<k && j<k &&i==j)
						h00[n*i+j]=1;
					else
						h00[n*i+j]=0;
				}
			}
		}
		double *h000=new double[n*n];
		h000=operatechengfa(Q,n,n,n,h00);
		for(int i=0;i<n;i++)
		{
			for(int j=0;j<n;j++)
			{
				Q[n*i+j]=h000[n*i+j];
			}
		}	
		int q=max(1,k);
		double *h1=new double[3*(n-q+1)];
		for(int i=0;i<3;i++)
		{
			for(int j=0;j<n-q+1;j++)
			{
				h1[(n-q+1)*i+j]=h[n*(k+i)+q-1+j];
			}
		}
		double *h2=new double[3*(n-q+1)];
		h2=operatechengfa(h0,3,3,n-q+1,h1);
		for(int i=k;i<k+3;i++)
		{
			for(int j=q-1;j<n;j++)
			{
				h[n*i+j]=h2[(n-q+1)*(i-k)+j-q+1];
			}
		}
		int r=min(k+4,n);
		double h3[3*r];
		for(int i=0;i<r;i++)
		{
			for(int j=k;j<k+3;j++)
			{
				h3[3*i+j-k]=h[n*i+j];
			}
		}
		double *h4=new double[3*r];
		h4=operatechengfa(h3,r,3,3,h0);
		for(int i=0;i<r;i++)
		{
			for(int j=k;j<k+3;j++)
			{
				h[n*i+j]=h4[3*i+j-k];
			}
		}
		x=h[n*(k+1)+k];
		y=h[n*(k+2)+k];
		if(k<n-3)
			z=h[n*(k+3)+k];
			delete []h0;
			delete []h1;
			delete []h2;
			delete []h4;
			delete []h00;
			delete []h000;
	}
	double xiangliang2[2]={x,y};
	double *h5=new double[4];
	h5=Householderjuzheng(xiangliang2,2);
	double *h7=new double[6];
	for(int i=0;i<2;i++)
	{
		for(int j=0;j<3;j++)
		{
			h7[3*i+j]=h[n*(i+n-2)+n-3+j];
		}
	}
	double *h6=new double[6];
	h6=operatechengfa(h5,2,2,3,h7);
	for(int i=0;i<2;i++)
	{
		for(int j=0;j<3;j++)
		{
			h[n*(n-2+i)+n-3+j]=h6[3*i+j];
		}
	}
	double *h8=new double[2*n];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<2;j++)
		{
			h8[2*i+j]=h[n*i+j+n-2];
		}
	}
	double *h9=new double[2*n];
	h9=operatechengfa(h8,n,2,2,h5);
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<2;j++)
		{
			h[n*i+j+n-2]=h9[2*i+j];
		}
	}
	double *h55=new double[n*n];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			if(i<n-2 && j<n-2 &&i==j)
				h55[n*i+j]=1;
			else if(i>=n-2 && j>=n-2)
				h55[n*i+j]=h5[2*(i-n+2)+j-n+2];
			else
				h55[n*i+j]=0;
		}
	}
	double *h555=new double[n*n];
	h555=operatechengfa(Q,n,n,n,h55);	
	return h555;
}
int select_m(double a[],int n)
{
	int m;
	int l;
	for(int i=n-1;i>=0;i--)
	{
		if(a[n*i+i-1]!=0)
		{
			if(a[n*(i-1)+i-2]!=0  && i>=1);
			{
				m=n-1-i;
				break;
			}
		}
	}
	return m;
}	
int select_l(double a[],int n)
{
	int m;
	int l=0;
	for(int i=n-1;i>=0;i--)
	{
		if(a[n*i+i-1]!=0)
		{
			if(a[n*(i-1)+i-2]!=0  && i>=1);
			{
				m=n-1-i;
				break;
			}
		}
	}
	for(int i=n-1-m;i>=0;i--)
	{
		if(a[n*i+i-1]==0 && i>0)
		{	
			l=i+1;
			break;
		}	
	}
	return l;
}
void yingshiQR(double h[],int n)
{
	for(int k=0;k<100;k++)
	{
		double *a=new double[n*n];
		double *b=new double[n*n];
		a=operatechengfa(h,n,n,n,h);
		for(int i=0;i<n;i++)
		{
			for(int j=0;j<n;j++)
			{
				if(i!=j)
					b[n*i+j]=a[n*i+j]-(h[n*n-1]+h[n*(n-2)+n-2])*h[n*i+j];
				else
					b[n*i+j]=a[n*i+j]-(h[n*n-1]+h[n*(n-2)+n-2])*h[n*i+j]+h[n*n-1]*h[n*(n-2)+n-2]-h[n*i+j]+h[n*n-2]*h[n*(n-2)+n-1];
			}
		}
		double *c=new double[n*n];
		c=QR_Q(b,n,n);
		double *d=new double[n*n];
		d=zhuanzhi(c,n,n);
		double *e=new double[n*n];
		e=operatechengfa(d,n,n,n,h);
		double *f=new double[n*n];
		f=operatechengfa(e,n,n,n,c);
		for(int i=0;i<n*n;i++)
		{
			h[i]=f[i];
		}
		delete []a;
		delete []b;
		delete []c;
		delete []d;
		delete []e;
		delete []f;
	}
}
double *duichengQR(double aa[],int n)
{
	double *a=new double[n*n];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			a[n*i+j]=aa[n*i+j];
		}
	}
	for(int k=0;k<n;k++)
	{
		double *b=new double[n*n];
		double *c=new double[n*n];
		b=QR_Q(a,n,n);
		c=QR_R(a,n,n);
		double *e=new double[n*n];
		e=operatechengfa(c,n,n,n,b);
		for(int i=0;i<n*n;i++)
		{
			a[i]=e[i];
		}
		delete []b;
		delete []c;
		delete []e;
	}
	return a;
}
double *duichengQR_Q(double aa[],int n)
{
	double a[n*n];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			a[n*i+j]=aa[n*i+j];
		}
	}
	double *h=new double[n*n];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			if(i==j)
				h[n*i+j]=1;
			else
				h[n*i+j]=0;
		}
	}
	for(int k=0;k<n;k++)
	{
		double *b=new double[n*n];
		double *c=new double[n*n];
		b=QR_Q(a,n,n);
		c=QR_R(a,n,n);
		double *e=new double[n*n];
		e=operatechengfa(c,n,n,n,b);
		for(int i=0;i<n*n;i++)
		{
			a[i]=e[i];
		}
		double *d=new double[n*n];
		d=operatechengfa(h,n,n,n,b);
		for(int i=0;i<n*n;i++)
		{
			h[i]=d[i];
		}
		delete []b;
		delete []c;
		delete []d;
		delete []e;
	}
	return h;
}
void yingshi1QR(double h[],int n)
{
	int m;
	int l;
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			if(abs1(h[n*i+j])<=(abs1(h[n*i+i])+abs1(h[n*j+j]))*pow(10,-8))
				h[n*i+j]=0;
		}
	}
	m=select_m(h,n);
	l=select_l(h,n);
	int k=0;
	while(m!=n && k<=10)
	{
		k++;
		double h12[(n-m-l)*l];
		for(int i=0;i<l;i++)
		{
			for(int j=0;j<n-m-l;j++)
			{
				h12[(n-m-l)*i+j]=h[n*i+j+l];
			}
		}
		double h22[(n-m-l)*(n-m-l)];
		for(int i=0;i<n-m-l;i++)
		{
			for(int j=0;j<n-m-l;j++)
			{
				h22[(n-m-l)*i+j]=h[n*(l+i)+j+l];
			}
		}
		double h23[(n-m-l)*m];
		for(int i=0;i<n-m-l;i++)
		{
			for(int j=0;j<m;j++)
			{
				h23[m*i+j]=h[n*(l+i)+j+n-m];
			}
		}
		double *p=new double[(n-l-m)*(n-l-m)];
		p=shuangchongbuQR_Q(h22,n-l-m);
		double *_h22=new double[(n-l-m)*(n-l-m)];
		_h22=shuangchongbuQR(h22,n-l-m);
		double *_h12=new double[(n-l-m)*l];
		_h12=operatechengfa(h12,l,n-m-l,n-l-m,p);
		double *pt=new double[(n-l-m)*(n-m-l)];
		pt=zhuanzhi(p,n-m-l,n-m-l);
		double *_h23=new double[(n-l-m)*m];
		_h23=operatechengfa(pt,n-m-l,n-m-l,m,h23);
		for(int i=0;i<n-m-l;i++)
		{
			for(int j=0;j<n-m-l;j++)
			{
				if(i>=l && j>=l && i<n-m && j<n-m)
					h[n*i+j]=_h22[(n-m-l)*(i-l)+j-l];
				else if(i<l && j>=l && j<n-m)
					h[n*i+j]=_h12[(n-m-l)*i+j-l];
				else if(i>=l && i<n-m && j>=n-m)
					h[n*i+j]=_h23[m*i+j-n+m];
			}
		}
		for(int i=0;i<n;i++)
		{
			for(int j=0;j<n;j++)
			{
				if(abs1(h[n*i+j])<=(abs1(h[n*i+i])+abs1(h[n*j+j]))*pow(10,-8))
					h[n*i+j]=0;
			}
		}		
		m=select_m(h,n);
		l=select_l(h,n);
		delete []p;
		delete []pt;
		delete []_h12;
		delete []_h22;
		delete []_h23;
	}
}
double *tezhenzhishibu(double a[],int n)
{
	double *aa=new double[n*n];
	for(int i=0;i<n*n;i++)
	{
		aa[i]=a[i];
	}
	yingshiQR(aa,n);
	double *b=new double[n];
	int k=0; 
	for(int i=0;i<n-1;i++)
	{
		if(abs1(aa[n*(i+1)+i])<=pow(10,-4))
		{
			b[k]=aa[n*i+i];
			k++;
		}
		else
		{
			b[k]=(aa[n*i+i]+aa[n*(i+1)+i+1])*0.5;
			b[k+1]=(aa[n*i+i]+aa[n*(i+1)+i+1])*0.5;
			k=k+2;
			i++;
		}
	}
	if(abs1(aa[n*n-2])<=pow(10,-4))
	{
		b[n-1]=aa[n*n-1];
	}
	return b;
}
double *tezhenzhixubu(double a[],int n)
{
	double *aa=new double[n*n];
	for(int i=0;i<n*n;i++)
	{
		aa[i]=a[i];
	}
	yingshiQR(aa,n);
	double *b=new double[n];
	int k=0;
	double d;
	double e;
	double *c=new double[n]; 
	for(int i=0;i<n-1;i++)
	{
		if(abs1(aa[n*(i+1)+i])<=pow(10,-4))
		{
			c[k]=0;
			k++;
		}
		else
		{
			d=aa[n*i+i]*aa[n*(i+1)+i+1]-aa[n*i+i+1]*aa[n*(i+1)+i];
			e=aa[n*i+i]+aa[n*(i+1)+i+1];
			c[k]=sqrt(abs1(4*d-e*e))*0.5;
			c[k+1]=-sqrt(abs1(4*d-e*e))*0.5;
			k=k+2;
			i++;
		}
	}
	if(abs1(aa[n*n-2])<=pow(10,-4))
	{
		c[n-1]=0;
	}
	return c;
}
double feiduijiao(double a[],int n)
{
	double s=0;
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			if(i!=j)
				s+=a[n*i+j]*a[n*i+j];
		}
	}
	return s;
}
int sgn1(double x)
{
	int s=1;
	if(x<0)
		s=-1;
	if(x>0)
		s=1;
	return s;
}
double *xunhuanJacobi(double aa[],int n)
{
	double *a=new double[n*n];
	for(int i=0;i<n*n;i++)
	{
		a[i]=aa[i];
	}
	double g=0.01;
	int k=0;
	double x;
	double t;
	double s;
	double c;
	int q;
	while(k<=1000000)
	{
		for(int i=0;i<n;i++)
		{
			for(int j=0;j<n;j++)
			{
				k++;
				if(abs1(a[n*i+j])>g && i!=j)
				{
					x=(a[n*j+j]-a[n*i+i])/a[n*i+j]*0.5;
					t=1.0*sgn1(x)/(abs1(x)+sqrt(1+x*x));
					c=1.0/sqrt(1+t*t);
					s=t*c;
					double b[2*n];
					for(int p=0;p<n;p++)
					{
						b[p]=c*a[n*p+i]-s*a[n*p+j];
						b[n+p]=s*a[n*p+i]+c*a[n*p+j];
					}
					b[j]=0;
					b[n+i]=0;
					b[i]=c*c*a[n*i+i]-2*s*c*a[n*i+j]+s*s*a[n*j+j];
					b[n+j]=s*s*a[n*i+i]+2*s*c*a[n*i+j]+c*c*a[n*j+j];
					for(int p=0;p<n;p++)
					{
						a[n*i+p]=b[p];
						a[n*j+p]=b[n+p];
					}
					for(int p=0;p<n;p++)
					{
						a[n*p+i]=a[n*i+p];
						a[n*p+j]=a[n*j+p];
					}
				}
			}
		}
		g=g/n;

	}
	return a;
}
double *xunhuanJacobijuzheng(double aa[],int n)
{
	double *a=new double[n*n];
	for(int i=0;i<n*n;i++)
	{
		a[i]=aa[i];
	}
	double *h=new double[n*n];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			if(i==j)
				h[n*i+j]=1;
			else
				h[n*i+j]=0;
		}
	}
	double g=0.01;
	int k=0;
	double x;
	double t;
	double s;
	double c;
	int q;
	while(k<=1000000)
	{
		for(int i=0;i<n;i++)
		{
			for(int j=i+1;j<n;j++)
			{
				k++;
				if(abs1(a[n*i+j])>g)
				{
					x=(a[n*j+j]-a[n*i+i])/a[n*i+j]*0.5;
					t=1.0*sgn1(x)/(abs1(x)+sqrt(1+x*x));
					c=1.0/sqrt(1+t*t);
					s=t*c;
					double b[2*n];
					for(int p=0;p<n;p++)
					{
						b[p]=c*a[n*p+i]-s*a[n*p+j];
						b[n+p]=s*a[n*p+i]+c*a[n*p+j];
					}
					b[j]=0;
					b[n+i]=0;
					b[i]=c*c*a[n*i+i]-2*s*c*a[n*i+j]+s*s*a[n*j+j];
					b[n+j]=s*s*a[n*i+i]+2*s*c*a[n*i+j]+c*c*a[n*j+j];
					for(int p=0;p<n;p++)
					{
						a[n*i+p]=b[p];
						a[n*j+p]=b[n+p];
					}
					for(int p=0;p<n;p++)
					{
						a[n*p+i]=a[n*i+p];
						a[n*p+j]=a[n*j+p];
					}
					for(int p=0;p<n;p++)
					{
						b[2*p]=c*h[n*p+i]-s*h[n*p+j];
						b[2*p+1]=s*h[n*p+i]+c*h[n*p+j];
					}
					for(int p=0;p<n;p++)
					{
						h[n*p+i]=b[2*p];
						h[n*p+j]=b[2*p+1];
					}
				}
			}
		}
		g=g/n;

	}
	return h;
}						

